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KEFOKT  SUMMARY 


This  Is  the  final  quarterly  technical  report  on  a one-year  contract 
to  extend  and  improve  an  existing  empirical  node 1 for  worldwide  behavior 
of  ionospherlcally  imposed  radio-wave  sclentlllation.  The  objectives  of 
the  project  were  (1)  to  improve  the  accuracy  of  model-based  calculations 
of  the  Intensity-scintillation  index  and  (2)  to  develop  a capability  for 
full  description  (from  the  point  of  view  of  engineering  applications)  of 
the  first-order,  complex-signal  statistics  that  characterize  the  trans- 
ionospheric  radio  communication  channel.  A follow-on  project  has  been 
initiated  to  extend  the  channel  model  to  include  second-order  signal 
statistics  in  the  temporal,  spatial,  and  spectral  domains.  The  first 
priority  in  the  follow-on  work  will  be  co  calculate  the  fluctuation 
spectre  of  relevant  signal  parameters.  In  addition,  it  is  Intended  to 
extend  validity  of  the  model  into  the  multiple-scatter  regime. 

The  stated  objectives  of  the  just-completed  project  were  accomplished 
by  starting  with  the  scintillation  model  developed  by  Fremouw  and  Rino 
(1973)  , adding  a geomagnetic-activity  dependence  for  the  behavior  of 
scintillation-producing  ionospheric  irregularities,  and  substituting  a 
more  accurate  and  general  scattering  theory  for  that  of  Briggs  and  Parkin 
(1963)  used  in  the  earlier  model.  The  results  of  these  efforts  have 
been  implemented  in  two  computer  programs.  The  first,  entitled  RFMOD, 
contains  t _ main  elements  of  the  scattering  theory,  the  morphological 
model  for  ionospheric  irregularity  strength  and  other  parameters,  and 
subroutines  for  calculating  the  geometry  and  other  relevant  quantities 
from  user-specif ied  inputs.  The  main  outputs  of  RFMOD  are  the  intensity- 
scintillation  index  S^,  which  is  the  fractional  nns  fluctuation  in  re- 
ceived signal  incensity  (Briggs  and  Parkin,  1963),  and  an  analogous  phase 
scintillation  index  $nas  (labled  PRMS) , which  is  the  rms  fluctuation  of 
phase . 
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Three  related  first-order,  signal-statistical  parameters  also  are 
output  from  RFMOD , to  be  used  as  inputs  to  the  second  program,  called 
DIST.  The  latter  program  permits  calculation  of  the  first-order  dis- 
tribution functions  of  amplitude  and  phase  for  selected  values  of 
and/or  $ras.  There  are  three  types  of  output  from  DIST:  values  of  the 
probability  density  functions  (pdf)  for  amplitude  and  for  phase,  and 
values  of  the  cumulative  probability  function  (cpf)  for  intensity  (or, 
equivalently,  for  amplitude  since  the  argument  is  output  in  dB  relative 
to  the  undisturbed  level).  The  cpf  output  constitutes  estimation  of 
fade  margin  required  to  combat  amplitude  scintillation. 

RFMOD  and  DIST  represent  interim  versions  of  evolving  transiono- 
spheric  communication-channel  codes.  As  stated  above,  capabilities 
will  be  added  for  calculating  second-order  statistics.  More  relevant 
in  the  near  term  is  that  there  has  been  no  direct  testing  of  calculations 
relating  to  phase  scintillation.  As  was  recognized  at  the  outset  of  the 
project,  data  are  not  available  at  present  that  are  directly  relevant 
to  testing  phase  statistics  unambiguously.  Useful  data  are  expected 
from  forthcoming  observations  of  multifrequency  coherent  satellite 
beacons,  and  such  tests  will  be  made.  The  forthcoming  data  may  also 
result  in  model  changes  that  will  further  alter  the  calculated  results 
for  amplitude  scintillation,  because  testing  to  date  has  still  necessarily 
relied  on  some  observations  reported  in  terms  of  a partly  subjective 
scintillation  index.  In  addition,  while  the  assumption  of  weak,  single 
scatter  contained  in  the  Briggs -and-Parkin  theory  has  been  partially 
removed,  a more  complete  correction  for  multiple  scatter  effects  has  yet 
to  be  incorporated. 

This  report  represents  a summary  of  the  work  performed  during  the 
past  year  and  is  inteded  as  a status  report  on  the  evolving  SRI/ARPA 
channel  model.  Section  II  describes  the  scattering  theory  employed  and 
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compares  some  of  Its  salient  points  with  those  of  other  theories.  That 
section  includes  some  material  reported  in  earlier  quarterly  reports, 
for  the  sake  of  presenting  a coherent  discussion  under  one  co'er,  but 
there  is  no  attempt  to  describe  the  theoretical  development  * " chrono- 
logical detail.  Section  III  contains  a brief  review  of  tl <*  known  mor- 
phology of  scintillation-producing  ionospheric  irregularities  as  modeled 
in  RFMOD.  A functional  description  of  the  programs  RFMOD  and  DIST  is 
presented  in  Section  IV,  along  with  instructions  for  their  use.  Finally, 
in  Section  V,  we  present  our  conclusions  on  the  utility  and  limitations 
of  the  channel  model  in  its  present  form  and  some  prospects  for  further 
improvement . 

II  SCATTERING  THEORY  EMPLOYED 

The  basic  approach  to  modeling  the  transionospheric  communication 
channel  is  to  treat  the  channel  as  a time-and-space-varying  linear  filter 
linking  the  transmitting  and  receiving  antennas  of  a communication  sys- 
tem. Slowly  varying  effects  (on  a time  scale  of  many  minutes  or  hours) , 
such  as  group  delay,  polarization  change,  and  dispersion  of  signals 
propagating  through  the  smooth  ionosphere,  can  be  described  by  well- 
known  deterministic  methods  and  are  not  included  in  this  work.  What  is 
treated  are  the  relatively  rapid  and  random  variations  (on  a time  scale 
of  fractions  of  a second  to  a minute  or  so)  in  signal  parameters  that 
arise  due  to  scattering  by  irregularities  in  ionospheric  electron  density 
and  that  are  referred  to  as  scintillations.  Because  the  channel  can  be 
treated  as  a linear  filter,  standard  Fourier  techniques  can  be  used  to 
apply  the  results  to  modulated  signals  if  a sufficiently  general  descrip 
tion  is  provided  of  the  signal  statistics  resulting  when  a continuous 
(CW)  is  passed  through  the  channel  (Fremouw,  1969) . In  this  section, 
we  discuss  the  first-order  statistics  of  this  basic  signal. 


Any  CW  signal  can  be  represented  by  a vector  on  the  complex  plane, 

such  as  the  phasor  E in  Figure  1.  In  this  figure,  E * ^ 4-  iE^j  is  meant 

to  represent  the  randomly  time-varying  complex  amplitude  of  the  signal 

output  from  an  antenna  terminal.  In  describing  the  first-order,  complex- 

signal  statistics  of  E,  we  are  assessing  the  probability  that  the  tip  of 

the  phasor  lies  within  an  elemental  area  on  the  complex  plain.  To  begin, 

we  write  E as  the  resultant  of  its  long-term  mean  (e)  and  a zero-mean, 

randomly  varying  component  E . We  choose  to  reference  all  pi ases  to 

s 

that  of  <e),  which  is  the  same  as  the  phase  that  E would  have  in  the 
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FIGURE  1 PHASOR  DIAGRAM  FOR  CONSIDERATION  OF  FIRST-ORDER 
SIGNAL  STATISTICS,  SHOWING  ELLIPTICAL  CONTOUR  OF 
EQUAL  PROBABILITY  PREDICTED  BY  THE  GAUSSIAN- 
STATISTICS  HYPOTHESIS 


absence  of  channel -imposed  perturbations  because  we  have  excluded  slow 
I channel  variations  from  the  model.  (Thus  (E  ) = 0,  and  (E  ) = T)  = 

j y x 

* (e).)  The  magnitude  of  (E)  (*  T0  is  less  than  the  amplitude  of  the  un- 

disturbed signal,  however,  because  some  of  the  transmitted  energy  has  been 

scattered  into  the  randomly  varying  component  E . 

s 
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To  Drocec^  with  the  analysis,  it  is  convenient  to  break  E into 

s 

its  own  real  and  imaginary  components,  S and  S , respectively.  The 

x y 

convenience  stems  in  part  from  the  fact  that  the  v<*rianc*-«  of  S and 

2 2 2 X 

S , a and  a , respectively,  and  their  covariance  a can  be  calculated 
y x y xy 

in  terms  of  ionosphe  r parameters.  These  are  useful  calculations  to 

2 2 

perform  because,  <it  > ..r  very  least,  the  sum  of  and  o represents  the 

power  contained  in  the  scattered  signal,  which  we  shall  denote  as  R . 

o 

That  is, 


R - <E  E*>  - 02  + a2  (1) 

o s s x y 

We  shall  normalize  all  our  calculations  to  unity  total  received  power, 

so  R will  represent  the  fraction  of  flux  scattered --i .e . , the  forward- 
o 

scattering  coefficient  of  the  ionosphere.  A complex  quantity  analogous 
to  the  real  that  it  is  useful  is  the  following  (Rino  and  Fremouw, 
1973a) : 


B = <E  E ) * (a2  - a2)  + 2i  a2  (2) 

o s s x y xy 

If  it  should  be  that  the  quadrature  components,  E and  E , are 

x y 

jointly  Gaussian  random  variables,  the  utility  of  the  three  variances 

defined  above--or  alternatively  of  R and  B --would  be  very  great.  Indeed, 

o o 

in  this  instance  the  first-order  statistics  of  the  received  signal  E are 

totally  defined  by  either  of  these  sets  of  three  parameters  (i.e.,  the 

three  real  variances  or  the  real  R^  and  the  real  and  imaginary  parts  of 

B ) . Most  scintillation  observations  have  been  only  in  terms  of  the 
o 

magnitude  of  E (i.e.,  amplitude  or  intensity  scintillation)  because  a 
sufficiently  stable  phase  reference  for  measuring  its  complex  value  has 
seldom  been  available  under  scintillation  conditions.  Thus,  there  has 
been  no  direct  determination  of  the  compl tx-signal  statistics.  (This 


data  deficiency  should  be  remedied  within  the  next  year  by  means  of 
careful  measurements  of  signals  from  the  VHF-UHF  coherent  beacons  on 
ATS-6  and  the  Transit  satellites  and,  under  a wider  range  of  ionospheric 
conditions,  by  means  of  observations  of  the  forthcoming  DNA-002  multi- 
frequency beacon.) 

The  Gaussian-statistics  hypothesis  (i.e.,  that  the  quadrature  com- 
ponents are  jointly  Gauss  an  variates)  has  been  tested  indirectly,  however, 
by  means  of  amplitude-scintillation  data.  Rino,  Livingston,  and  Whitney 
(1975)  assessed  this  hypothesis  and  the  competing  one  of  log-normal 
statistics  (according  to  which  the  phase  and  the  logarithm  of  amplitude 
are  jointly  Gaussian  variates)  by  performing  chi-square,  goodness -of-f it 
tests  of  amplitude  histograms  to  amplitude  pdf's  calculated  from  the  two 
hypotheses.  Good  fits  were  found  in  general  for  both  hypotheses,  which 
tend  toward  identical  pdf's  as  the  scintillation  index  decreases,  but 
the  Gaussian  hypothesis  yielded  better  fits  to  11  out  of  the  12  data 
sets  tested.  This  result  is  fortunate  because  it  permits  presentation 
of  scintillation  calculations  in  a form  that  will  be  rather  familiar  to 
communications  engineers.  In  the  remainder  of  this  report,  we  proceed 
on  the  basis  of  the  Gaussian  signal-statistics  hypothesis. 

For  Gaussian  signal  statistics,  a contour  of  equal  probability  for 

the  tip  of  the  E-phasor  on  the  complex  plane  is  an  ellipse  (Beckman  and 

Splzzichlno,  1963),  such  as  the  one  in  Figure  1.  The  analytical  task  of 

first-order,  signal-statistical  channel  modeling,  then,  is  to  relate  the 

parameters  of  such  an  ellipse  to  parameters  of  the  channel  (scattering 

irregularities  in  the  ionosphere  in  the  present  case)  . The  probability 

ellipse  is  characterized  by  its  size,  eccentricity,  and  orientation  on  the 

complex  plane.  Clearly,  R (being  the  fraction  of  flux  scattered)  dic- 

o 

tates  the  size  of  the  ellipse.  The  eccentricity  is  dictated  jointly  by 

T.  and  the  magnitude  of  B , and  the  orientation  is  controlled  by  the  phase 
o o 

angle  of  B^.  Specifically,  the  axial  ratio  of  the  probability  ellipse  is 
given  by 


6 


n 1/2 


r 


(3) 


and  its  orientation  angle,  6 (shown  on  Figure  1),  is  given  by 


6 


(4) 


It  will  be  recognized  from  Eqs.  (1)  through  (4)  and  the  accompanying 

discussion  that  there  are  several  sets  of  three  parameters  that  can  be 

used  to  describe  the  first-order,  complex-signal  statistics.  For  reasons 

of  convenience,  which  will  became  clear  in  the  ensuing  discussion,  we 

choose  to  work  with  the  parameters  R , |b  |/R  , and  6.  As  an  example  "f 

o o o 

the  relationship  of  these  three  parameters  to  a quantity  of  direct 
applications  interest,  we  have  plotted  the  intensity-scintillation  index 
in  Figure  2 from  the  following  formula  derived  by  Nakagami  (1960) : 


S - 2R  (1  - R ) 1 • 

4 o o 


( 


cos  26 


)■<{'•¥ 


(5) 


We  note  that  in  Briggs'  and  Parkin's  theory,  S,  is  calculated 

4 

approximately  as 


I 


? 

t 

f 

I 


(6) 
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FIGURE  2 SCINTILLATION  INDE'i  S4  AS  A FUNCTION  OF  ELONGATION  AND 
ORIENTATION  OF  THE  EQUAL-PROBABILITY  ELLIPSE,  FOR  FOUR 
VALUES  OF  IONOSPHEW.C  SCATTERING  COEFFICIENT,  Rq 


The  following  three  relationships  arise  from  F.qs . (1)  and  (2): 


2 

a 

xy 


8 in  (26) 


(9) 


Thus,  Eq.  (6)  contains  only  the  contribution  of  the  in-phase  scattered 
component  to  S^.  It  arises  from  equating  fluctuations  in  to 

fluctuations  in  the  real  amplitude  A and  is  accurate  only  for  very  weak 
scatter  [i.e.,  when  the  higher-order  terms  in  Eq.  (5)  are  negligible]. 

A significant  improvement  in  the  present  scintillation  model  is  retention 
of  all  terms  in  Eq.  (5) , which  accounts  for  the  contribution  of  fluctua- 
tions in  the  phase-quadrature  component  E^(“S  ) to  amplitude  scintilla- 
tion (as  well  as  for  the  influence  of  the  correlation  between  S and  S ) . 

x y 

Other  differences  between  the  old  scattering  theory  and  the  new 

arise  in  the  manner  of  calculating  R and  B . The  basic  approach  used 

o o 

in  the  present  work  is  that  described  by  Rino  and  Fremouw  (1973b) , in 
which  the  scattering  geometry  is  treated  quite  rigorously.  Briggs  and 
Parkin  simplified  the  geometry  by  assuming  that  the  ionospheric  scattering 
layer  effectively  can  be  tilted  from  the  horizontal  so  that  its  boundaries 
are  always  normal  tc  the  incident  radio-wave  propagation  vector.  This 
approximation  breaks  down  for  large  zenith  angles  but  it  avoids  a great 
deal  of  complication  in  derivation  of  the  scattering  equations.  Some  of 
this  complication  was  avoided  by  Rino  and  Fremouw  (1973b)  by  assuming 
zero  correlation  in  electron-density  fluctuation  along  the  vertical,  which 
has  negligible  effect  on  the  scattering  equations  for  nearly  vertical 
incidence.  Again,  however,  the  effect  of  the  assumption  becomes  un- 
acceptable as  the  incidence  angle  increases. 

During  the  last  quarter  of  this  contract,  the  theory  was  recast  to 
include  the  scattering  effect  of  structure  in  all  three  dimensions  of 
the  ionosphere.  An  outline  of  the  derivation  is  given  here,  together 
with  the  resulting  equations  for  describing  first-order  signal  statistics, 


9 


which  are  valid  over  a greater  range  of  incidence  angles  than  their 
counterparts  in  Quarterly  Technical  Reports  1 and  2. 

The  scattering  coefficient  R represents  a special  case  o-f  the 

o 

spatial  autocorrelation  function  of  the  scattered  signal,  which  is 
defined  as 


R(P)  m <E  (r)  E*  (r  + p)  > ; (10) 

s s 

namely,  the  case  in  which  the  spatial  lag  parameter  P is  zero.  A similar 

definition,  calculated  without  taking  the  complex  conjugate,  can  be  made 

of  a quantity  B(p)  , which  reduces  B in  the  same  special  case  (l.e., 

° 

zero  lag) . The  approach  to  evaluating  R(P)  and  B(p)  is  to  express  E 

s 

as  the  scattered  field  calculated  by  means  of  the  first  Born  approxi- 
mation, with  the  scattering  medium  represented  by  a spatial  spectrum  of 
ionospheric  structure.  The  field  is  calculated  by  inte  ating  the 
contribution  from  different  heights  within  the  scattering  layer,  con- 
sidering the  diffraction  effects  arising  in  propagation  to  the  observing 
plane  but  ignoring  effects  of  multiple  scatter.  When  the  complex  field 
is  multiplied  by  its  complex  conjugate,  the  effects  cancel. 

The  result  is  that  R(p)  can  be  expressed  as  the  two-dimensional 
Fourier  transform  of  a height-smeared  version  of  the  ionosphere's  spatial 
spectrum  and  B(p)  can  be  expressed  as  the  transform  of  that  spectrum 
times  a propagation  factor.  Setting  the  lag  parameter  to  zero  is  tanta- 
mount to  integrating  over  the  spectra  in  the  observing  plane  (as  opposed 

to  calculating  the  Fourier  transforms) , with  R being  an  integral  over 

o 

the  scattered  field's  angular  power  spectrum  and  B being  an  integral 

o 

over  a complex  spectrum  that  incorporates  diffraction  effects. 

Probably  more  important  for  first-order  statistics  than  the  approach 

to  deriving  R and  B is  the  form  of  the  spatial  spectrum  used  in 
o o 
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their  calculation.  In  evaluating  R and  the  second  factor  in  Eq.  (6), 

o 

Briggs  and  Parkin  assumed  a Gaussian -shaped  spatial  spectrum  for  the 
ionosphere.  While  convenient  because  of  its  Fourier- transformation 
characteristics,  a Gaussian  spectrum  does  not  appear  to  represent 
ionospheric  structure  very  realistically  (Rufenach,  1972;  Dyson  et  al . , 
1974)  . 

A better  description  is  afforded  by  a power  law  with  a low-frequency 

cutoff.  Whan  normalized  so  that  its  three-dimensional  integral  is  unity 

2 

times  the  spatial  variance  ((AN)  ) of  electron  density,  such  a spectrum 

3 2 

[units  of  (electrons/m  ) per  unit  three-dimensional  spatial  frequency, 

OS 

if]  may  be  written  as 


K(tC)  - 8TT 


3/2  + 1/2) 

r<v  - i) 


3 

or  a 


<(AN)  ) 


1 + 


[<*sO)]2 


v + 1/2 


(ID 


where  v is  the  spectral  index  defining  the  power  law,  a is  scale-size 
parameter  controlling  the  low-frequency  cutoff,  a is  the  ratio  of 
Irregularity  scale  size  along  the  geomagnetic  field  to  that  across  the 
field,  and  the  three  dimensional  shape  and  orientation  of  the  irregulari- 
ties (Budden,  1963)  , are  described  by 


S2(*>  - A(t)  K2  + K2  9-  B(t)*2  - 2C(*)K  K 
x y z x z 


(12) 


where 


A(t) 


lL±  92ih 

BO) 


03) 


B(  ♦) 


2 2 2 
cos  y + a sin  t 


(14) 
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2 

C(  *)  * (1  - a ) cos  i sin  if 


(15) 


and  t * magnetic  dip. 


The  coordinate  system  used  is  such  that  is  directed  to  the  geomagnetic 
north,  * to  the  geomagnetic  east,  and  K downward. 

y z 

Using  the  above  spectrum,  one  obtains  the  following  expression  for 

R : 

o 


R 

o 


. 2 ,2  T(v  ■ 1/2)  a a sec  6 _ , 2, 

V"  r.  x r(v  r -7=5-  L«iN>  > 


(16) 


where 


and 


r " classical  electron  radius, 

e 

X ■ oboerving  wavelength, 

L “ thickness  of  scattering  layer, 

A'  « 1/2[b'  + D], 

C'  « 1/21.D ' - D], 

D'  - 1 + B(*  + 90°)  + 2C(+)  tan  9 cos  0 + B(f)  tin2  9, 

D - /Id'  - 2(1  + sin^  0)]2  + (B')2  , 

b'  - B(\|i)  tan  9 sin  20  + 2C(>li)  tan  9 sin  0, 

a 

6 « incidence  angle  on  scattering  Icyer, 

0 * propagation  azimuth  relative  to  geomagnetic  north  > *■ 

height  of  scattering  layer. 


Eq.  (16)  was  obtained  by  integrating  analytically  over  the  angular 

power  spectrum  of  the  scattered  wave  arriving  at  the  observing  plane, 

as  described  above.  The  integral  over  the  complex  spectrum  associated 

with  B cannot  be  performed  analytically.  However,  the  real  and  imaginary 
o 

integrals  involved  have  been  solved  numerically  for  a large  range  of 


12 


geometrical  situations  described  by  the  parameters  a,  \|t , 6,  and  Study 
of  the  numerical  results  has  led  to  approximations  suitable  for  implemen- 
tation in  an  operating  computer  code  such  as  RFMOD. 

Simplifications  have  been  found  by  converting  the  real  and  imaginary 

integrals  to  the  magnitude  and  phase  angle  of  B^  and  then  to  the  quantities 

|b  1 /R  and  6 [the  latter  via  Eq.  (4)].  The  simplification  arises  from 
o o 

an  empirically  noted  approximate  relationship  between  the  latter  two 

quantities t as  described  in  detail  in  Quarterly  Technical  Reports  1 and  2. 

It  permits  approximate  description  of  the  complete  first-order,  complex- 

signal  statistics  from  solution  of  the  real  integral  alone.  What  is 

2 

actually  calculated  is  o /R  , which  is  related  to  Re[B  } by  Eq.  (7). 

x o o 

Two  approximations  are  involved  in  formulating  this  part  of  the 

problem.  The  first- -neglect  of  a diffraction  effect  related  to  layer 

thickness--appears  numerically  to  be  valid  essentially  anywhere  outside 

the  scattering  layer.  The  second  holds  so  long  as  the  spectrum  of 

irregularities  contributing  to  fluctuations  is  cut  off  at  the  low- 

frequency  end  by  diffraction  effects  rather  than  by  an  inherent  cutoff 

of  the  ionospheric  spatial  spectrum.  Both  in-situ  (Dyson  et  al . , i974) 

and  scintillation  observations  (Rufenach,  1972;  Rino,  Livingston,  and 

Whitney,  1975)  support  this  view  for  all  frequencies  of  interest.  We 

have,  therefore,  proceeded  on  the  basis  of  this  "near-zone"  approximation, 

but  we  have  incorporated  a test  in  RFMOD  to  safeguard  against  violation 

of  the  assumption  underlying  it.  On  the  basis  of  these  two  approximations, 
2 

one  may  write  a /R  as  follovs: 
x o 


I'-r 


sin2  [q2  7(6.0  ')  J 


dtp  dq 


07) 


UP')  7 coa2  tp  + sin2  cp]V  + 1/2 


where  \/\Z * Fresnel-zone  radius  at  scattering  layer 

7(0,0')  - [l  + tan  6 cos2  (0'  - 9)3^  , 


and 


(18) 


8'  - v/a '/c ' 

0 ' * 0 - ” tan 


B 


d'  - 2(1  + sin  0) 


(19) 


The  main  diffraction  effect  dictating  the  first-order  statistics  of 
a scattered  signal  is  controlled  by  the  Fresnel-distance  parameter 
appearing  in  parentheses  in  front  of  the  integral  in  Eq.  (17)  . The 
integral  Itself  contains  only  geometric  factors.  Although  the  integral 
canno.t  be  solved  analytically,  it  has  been  possible  to  find  a useful 
algebraic  approximation  to  it.  The  approach  was  to  numerically  evaluate 
the  integral,  which  was  named  F(0',0,  0,),  for  selected  points  in  a 
three-space  of  0',  0,  and  0',  and  then  to  fit  judiciously  chosen 
algebraic  functions  to  the  points.  The  resulting  functions  then%ere 
incorporated  into  RFMOD.  The  procedure  was  as  follows. 

First,  the  integral  was  solved  numerically  for  0 ' * 1,  2,  5,  10,  and 
20,  with  6 ranging  from  0°  through  80°  (which  exce.'dt  the  maximum  in- 
cidence angle  of  a trans ionospheric  communicatiot.  \ Vi  on  the  F layer, 
by  virture  of  the  ionosphere's  curvature),  and  with  t stepped  between 
0°  and  90°  (since  the  0 dependence  is  symmetrical  about  -0°;  . After 
the  results  were  plotted,  inspection  of  the  curves  suggested  behavior 
of  the  following  form: 

F(0  ' , 6,  0')  * ja  (0',  6)  + b(0',  9)  + (b  - a)  cos  20',  (20) 
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where  a ■ the  value  of  the  Integral  when  0'  «*  tt/2,  and  b ■ the  value 
when  0 ' * 0 . 

The  values  of  a and  b were  then  plotted  as  functions  of  6 for  given 
values  of  P-.  Again  the  plots  suggested  analytic  fits  of  the  form 

a(p',G)  = C(P')  8ec  0m^  * 1 (21) 

b(P',9)  - c(P')  sec  9n(P  ' (22) 

where  c(P  ')  * the  value  of  the  function  when  Q ■ 0 / - 0 and  where  the 
exponents  m(P  ')  and  n(P ')  were  chosen  by  forcing  Eqs.  (21)  and  (22)  to 
match  the  numerical  value  of  the  integral  when  9 * 70°,  which  is  about 
its  largest  value  for  penetration  of  the  F layer.  The  values  of  c,  m, 
and  n were  then  employed  in  a one -dimensional , least-squares  fitting 
routine  to  obtain  best-fit  polynomial  expansions  of  c(0  ') , ra(P  ') , and 
n(3  ')  . 

The  results  of  the  above  procedure  were  combined  into  an  algebraic 
approximation  to  the  integral  in  Eq.  (17) , and  the  approximation  was 
then  used  to  calculate  curves  of  FO',  9,  0 ')  over  the  ranges  of  9 and 
0'  and  for  the  values  of  P'  originally  employed  in  the  numerical  in- 
tegration. Finally,  similar  curves  were  calculated  for  other  (both 
Intermediate  and  more  extreme)  values  of  and  numerical  integrations 
were  performed  for  these  new  cases.  The  entire  collections  of  numerical- 
integration  and  algebraic-approximation  results  were  then  compared  to 
determine  the  efficacy  of  the  approximate  expression.  An  example  of  the 
comparison  is  presented  in  Figure  3 for  a value  of  P ‘ not  involved  in 
development  of  the  approximation;  the  smooth  curves  represent  the 
algebraic  fit  and  the  individual  points  are  results  of  the  numerical 
integration.  The  fit  is  considered  quite  adequate  for  application  to 
scintillation  modeling.  Thus,  RFMOD  contains  the  following  expression: 


» 

; 
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X 


FO  6,  0;)  , 


(23) 


where  FO*,  9,  0 ')  is  calculated  from  Eq.  (20). 


LA-341S-S 


FIGURE  3 COMPARISON  OF  ALGEBRAIC  APPROXIMATION  (Eq.  20)  TO  VALUES  OF  THE 
INTEGRAL  IN  EQ.  (17) 


It  is  clear  from  Eqs . (7),  (8),  and  (9),  that  the  first-order 

statistical  problem  would  be  completely  solved  at  this  point  if  one 

2 

could  find  an  approximation  to  o in  some  manner  analogous  to  that 
2 2 Xy  2 

above  for  a (a  being  simply  R - o ) . However,  no  such  approximate 


solution  to  the  imaginary  part  of  the  B integral  has  been  found,  and 

o 

numerical  integration  each  time  RFMOD  is  run  seems  unacceptably  inefficient. 
Fortunately,  as  mentioned  above,  synthesis  of  the  numerical  integration 
results  that  have  been  obtained  has  pointed  to  another  simple  procedure. 
Specifically,  when  Ib^I/R^  and  6 were  formed  from  the  real  and  imaginary 
integrals  for  a number  of  geometrical  situations,  it  was  found  that  the 
relationship  between  |b  |/F^  and  6 was  not  very  dependent  on  geometry 
and  was  well  approximated  by  the  following  simple  formula: 


6 - j(l  - |boI/Ro).  (24) 

Not  only  is  the  relationship  between  6 and  |b  | /R  rather  insensitive 

o o 

to  geometry,  but  also  the  first-order  signal  statistics  are  not  very 

sensitive  to  that  relationship.  For  first-order  statistics,  the  main 

effect  of  diffraction  is  a simultaneous  rotation  (increasing  6)  and 

circularization  (decreasing  |b  | /R  ) of  the  correlation  ellipse  illus- 

o o 

trated  in  Figure  1,  as  the  scattered  wave  propagates  with  increasing  Z in 

Eq.  (23).  What  is  important  is  the  rate  at  which  this  approach  to  Rice 

statistics  (|b  | /R  * 0 and  6 = 45°)  takes  place  rather  than  the  precise 
o o 

relation  between  ellipse  elongation  and  eccentricity.  Thus,  we  have 

proceeded  on  the  basis  of  Eq.  (24),  which,  when  combined  with  Eq.  (7), 

2 

yields  the  following  transcendental  equation  relating  Ib  | /R  to  o /R  : 

o o x o 


Eq.  (25)  is  fit  very  well,  especially  in  the  near  and  intermediate 
zones,  by  the  following  approximation: 


( 26) 
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1 


I 
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Inserting  Eqs.  (7)  and  (26)  into  Eq.  (5),  we  obtain  the  following  very 
useful  formula  for  scintillation  index: 


While  containing  some  approximations  in  its  third  term,  Eq.  (27)  is  a 

considerable  improvement  over  the  formula  of  Briggs  and  Parkin  for 

[Eq.  (6)],  which  totally  neglected  both  the  second  and  third  terms.  The 

improvement  increases  as  the  scattering  becomes  more  intense  (increasing 

R ). 

o 

One  reason  for  separating  the  first  two  terms  in  Eq.  (27)  into  two 

factors  is  that  R is  amenable  to  a simple  multiple-scatter  correction, 
o 

The  theory  discussed  so  far  is  limited  to  the  regime  of  weak,  single 
scatter  for  the  entire  scattering  region  (i.e.,  Rq  « 1).  However,  an 
important  effect  of  multiple  scatter,  namely  extinction  of  the  non- 
scattered  signal  component  (represented  by  (e)  in  Figure  1)  , can  be 
accounted  for  simply  by  requiring  conservation  of  total  power  flux  while 
applying  the  first  Born  approximation  to  elemental  layers  within  the 
scattering  region.  From  the  work  of  Fejer  (1953) , Bramley  (1955)  , 
Fremouw  (1968),  and  Rino  (1975),  it  follows  that  extinction  can  be 
accounted  for  by  replacing  the  single-scattering  coefficient  calculated 
from  Eq.  (16)  by  a multiple-scattering  coefficient  defined  by 


R = 1 
o 


-Rr 


(28) 


It  is  R/  that  is  used  (in  place  of  R ) in  Eq.  (27)  for  calculation 
o o 

of  scintillation  index,  which  is  output  as  S4,  in  RFMOD;  this  multiple- 
scattering coefficient  (R ')  is  itself  output  from  RFMOD  as  RZRO.  An 

o 
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important  task  £or  the  second  contract  year  is  to  find  an  analogous 

2 

multiple-scatter  correction  for  o /R  and  the  related  quantities  IB  /R 

x o o o 

and  6.  Presently,  the  latter  two  parameters  are  calculated  from  Eqs. 

(26)  and  (24),  respectively,  on  the  basis  of  Eq.  (23).  They  are  output 

from  RFMOD  as  BMAGR  and  DELTA,  respectively,  for  use  as  inputs  to  the 

second  channel-model  code,  DIST.  A single-scatter  approximation  warning 

is  printed  if  R a 1. 

o 

Once  R ' , |b  | /R  , and  6 are  calculated,  then  application  of  the 
o o o 

Gaussian-statistics  hypothesis  r nits  calculation  of  any  desired  first- 
order,  signal-statistical  parameter.  The  calculation  of  the  intensity- 
scintillation  index, 


s A «A2  - <*V>1/2 

u2> 


(29) 


as  described  above,  is  an  example.  An  analogous  phase-scintillation 
index  may  be  defined  as 


i _ ((«  - <*>)2>1/2 

rms 


- (<$2>  - <$>V/2 


(30) 

(31) 


Given  the  phase  reference  chosen  in  Figure  1,  it  is  well  known  (e.g. 

Beckman  and  Spizzichino,  1963)  that  under  the  postulate  of  Gaussian 

statistics,  the  joint  probability  density  function  (pdf)  for  the  real, 

E , and  imaginary,  E , parts  of  the  received  signal's  complex  amplitude, 
x y 


E, 


is  giver  by 


P 


(E  ,E  ) - 


xy  x y 


1 


„ 2 If  2 

2tto  a VI  - p 
x y 


exp 


2p 


(E  - T])E  E 

x y 

a a 
x y 


2 1 


y J 
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(32) 


where  p » o2  /o  o end  where  T|  mey  be  celculeted  from  T]  * >/l  - R *nd 
xy  x y ° 

the  o's  from  Eqe.  (7),  (8),  end  (9).  The  indlviduel  pdf'a  for  amplitude 
end  for  phese  cen  be  calculated  by  changing  verieblea  in  Eq.  (32)  fr<xn 
c exponent  to  poler  fora  end  then  integreting  over  the  other  verieble. 

Thet  ie, 

2tt 

p(A)  - £ A Pxy(A  cos  *»  A sin  d*  <33) 

m 

and  P(*)  ■[  A p (A  cos  i,  A sin  *)  dA  (3^) 

•k  xy 

Analytic  solutions  to  the  integrals  in  Eqs.  (33)  and  (34)  have  been 

known  for  some  time  in  special  cases  (e.g.,  for  Rice  statistics,  in  which 

a - 0 and  a - a ) . but  not  so  for  the  generalised  Gaussian  statistics 
xy  x y 

needed  to  describe  scintillation.  Recently,  however,  Hatfield  (1975) 
has  found  the  following  general  analytic  solution  for  Eq.  (34): 


exp 


p(4)  - 


o 

X 


*-2(i  - g*) 


2ttct  a >/l  - p2  f}(i) 
x y 1 


(1 


-p2) 


exp 


- f2(«) 


2(1  - P ) 


+ f 2(«)  2TT^2'  P~^~  erfc 


- f2(«) 

y/m  - P2) 


(35) 


where  f (t)  - \/cos2  i/a2  - 2p  cos  i sin  i/(o  a ) + sin2  i/a 


x y 


(36) 


and  f„(*)  - 'H  [ecs  i/a2  - p sin  t/(o  o)]/f  (*) 
2 X X x y x 


(37) 
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Unfortunately,  there  is  no  known  analytic  solution  for  the  moments  of 
the  distribution  in  Eq.  (35).  Therefore,  the  distribution  itself  has 
been  incorporated  in  RFMOD,  together  with  a numerical  Integration  routine 
for  calculating  its  first  and  second  moments,  as  follows: 


ft  p («)  dt 


<$2)  - f t2  p ( t ) dft 
*o 


The  results  of  Eqs.  (38)  and  (39)  are  then  inserted  into  Eq.  (31)  for 
calculation  of  the  phase-scintillation  index,  which  is  output  as  PRMS. 

As  described  above,  the  signal-statistical  outputs  of  RFMOD  are 

R ' , Ib  I /R  , 6,  S , and  t . The  last  two  outputs  constitute  an  intensity- 
o o o 4 rms 

scintillation  index  and  a phase-scintillation  index,  respectively,  and 
represent  our  judgment  as  to  the  first-order  statistical  parameters  most 
likely  to  be  sought  by  a user  for  synoptic  purposes.  The  other  three 
outputs  are  provided  for  use  in  posing  more  detailed  questions  about 
first-order  signal  statistics  in  particular  cases,  by  means  of  the 
DIST  code . 

The  DIST  code  is  essentially  a device  for  calculating  and  outputing 
the  amplitude  and  phase  pdf's  from  inputs  of  R ' , |b  | /R  , and  6.  For 

COO 

phase,  it  is  simply  a matter  of  evaluating  Eq.  (35)  for  the  given  input 
parameters  and  outputing  the  result  as  P(PHI) . Unfortunately,  there  is 
no  known  analytical  solution  to  Eq.  (33)  analogous  to  Hatfield's  phase 
distribution.  Thus,  it  is  necessary  to  perform  the  phase  integration 
over  the  joint  distribution  (32)  numerically.  This  integration  is  per- 
formed in  DIST,  and  the  result  is  output  as  P(A)  versus  A (the  latter  in 


unit*  of  undisturbed  amplitude).  In  addition,,  the  following  integration 
also  is  performed  numerically: 

V 

P(A  < V)  - f P(A)  dA  . (40) 

0 

This  resulting  cumulative  probability  function  (cpf)  of  amplitude  is 
then  output  as  CPF  versus  20  log  V (labeled  dB) . Since  the  independent 
variable  is  output  in  decibels,  CPF  may  be  thought  of  as  the  cumulative 
distribution  of  intensity  (1  - A^)  as  well  as  of  amplitude  (A) ; it 
represents  the  link  margin  needed  to  mitigate  scintillation  fading  to  a 
desired  reliability  threshold. 


) 


1 

I 


Ill  SCINTILLATION  MORPHOLOGY  MODELED 


The  major  work  performed  on  thin  contract  was  development  and 
Implementation  of  the  generalized  scattering  theory  described  in  Section 
II.  The  other  important  aspect  of  modeling  the  trans ionospheric  scintil- 
lation channel  is  description  of  the  ionosphere's  scattering  structure  as 
a function  of  several  variables.  The  structural  characteristics  that 
must  be  d%scriDed  pertain  to  the  strength,  size,  and  shape  of  the  radio- 
wave  -scattering  irregularities  and  the  thickness  of  the  scattering  layer. 
All  the  relevant  parameters  appear  in  Eq.  (16).  The  shape  of  the  irregu- 
larities is  described  by  the  ionospheric  spatial  spectrum  [Eas.  (11)  and 
(12)],  which  is  characterized  mainly  by  the  power-law  spectral  index,  v, 
and  by  the  ratio  a,  of  scale-sisse  pa^am^ter  along  the  geomagnetic  field 
to  that  across  It.  The  size  of  the  irregularities  is  scaled  by  the 
cross-field  scale-size  parameter,  a.  The  strength  of  the  irregularities 
is  described  by  the  spatial  variance  of  electron  density  <(£N)^),  and 
the  thickness  of  the  scattering  region  is  L. 

In-situ  measurements  of  the  ionospheric  spatial  spectrum  (Dyson  et 
al.,  1974)  show  a'remarkably  consistent  spectral  index,  corresponding 
very  nearly  to  the  Kolmorov  value  of  4/3  for  v,  under  most  ionospheric 
conditions;  we  have  adopted  this  value  as  a constant  in  the  model.  The 
irregularity  axial  ratio,  a,  also  has  been  taken  as  constant,  with  a 
value  of  10,  as  a result  of  the  data  review  performed  by  Fremouw  and 
Rino  (1973) . There  is  no  direct  measurement  of  the  outer  scale,  a;  it  is 
simply  known  to  be  large  compared  with  a km  or  so  (Rufenach,  1972;  Dyson 
.t  al,,  1974;  Sagalyn  et  al . , 1974).  We  have  found  that  a * 3 km  pro- 
duces a good  fit  to  the  amplitude-scintillation  data  analyzed  by  Rino, 
Livingston,  and  Whitney  (1975)  , and  this  value  is  currently  employed  as 
a constant  in  the  model. 
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The  remaining  two  ionospheric  parameters,  ((AN) 2)  and  L,  appear  in 
the  scintillation  calculations  only  as  a product.  We  have  set  L *■  100  km 
as  a representative  value  and  have  expressed  the  behavior  of  scintillation 
in  terms  of  the  following  empirical  model  for  the  rms  electron -density 
variation,  <(AN)2>1/2: 


W>1/2 


AN  a,  T ,D,R)  + AN  (X,T)  + AN  (X.T.K  ) 
e m n p 


+ AN  (X,T,K  ) el/m^ 
* P 


where  X 
T 
D 
R 
K 

P 


geomagnetic  latitude  in  degrees 
time  of  day  in  hours, 
day  of  the  year  out  of  365, 
sunspot  number, 
planetary  magnetic  index. 


(41) 


The  terms  in  Eq.  (41)  are 


AN  - (1.2  X 10  )(1  + 0.05R) 


1 - 0.25  cos  tt 


+ exp 


3 if  T 4 22.5 
D 16  if  T > 22.5  (42) 


AN  - (4.4  X io8)(i  + o.: 


m 


TTT  \ 

/x  - 32.5V  1 

COS  12/ 

exp 

A 15  / 

1 

L \ / J 

AN  - (2.0  X 10  ) 
n 


1 + ^rf 


X - 68.5  + 1.82  Kp  + (3  + 0.55  cos  (nT/‘12,J 


18.7  - 0.36  Kp  - (1  + 0.11  Kp)  cos  (ttT/12) 


(43) 


(44) 


(45) 


«nd  AN  - (1.1  X 10  ) 


K exp 
P 


X - 77  + 1.82  Kp  +4  cos  (TrT/12y 


0.32  Kp 


I1 


Eqs.  (42)  and  (43)  describe,  respectively,  the  equatorial  and  mid- 
latitude  behavior  of  scintillation,  and  the  high-latitude  behavior  is 
described  by  Eqs.  (44)  and  (45).  All  four  terms  of  Eq.  (41)  contain  the 
we 11 -established  nighttime  maximum  of  scintillation.  The  equatorial  term 
also  includes  a seasonal  variation,  displaying  equinoctial  maxima.  The 
most  complicated  behavior  occurs  at  high  latitudes,  above  the  "scintil- 
lation boundary"  (Aarons  et  al.,  1969).  Equation  (44)  essentially 
describes  the  behavior  of  the  boundary  location  as  a function  of  time  of 
day  and  geomagnetic  activity,  and  Eq.  (45)  accounts  for  scintillation 
that  is  directly  related  to  auroral  activity. 

In  the  Fremouv-and-Rino  1973  model,  the  relationship  of  scintilla- 
tion to  geomagnetic  activity  was  not  Included.  The  only  solar- terrestrial 
observable  employed  In  the  earlier  model  was  the  sunspot  number,  which 
was  Intended  as  a measure  of  epoch  during  a solar  cycle.  However,  it  is 
clear  that  the  scintilla* ton  boundary  migrates  equatorward,  along  with 
the  auroral  oval,  with  increasing  planetary  magnetic  activity  (Aarons 
et  al.,  1963;  Evans,  1974).  This  was  known  at  the  time  the  1973  model 
was  developed,  but  the  dependence  was  excluded  to  avoid  double  accounting 
for  other  aspects  of  scintillation  behavior,  such  as  seasonal  and  solar- 
cycle  variation. 

Pope  (1974)  modified  our  high-latitude  term  to  include  a K de- 

P 

pendence  without  checking  on  the  possibility  of  seasonal  double  accounting. 
We  have  had  five  years  of  values  plotted  and  have  found  no  troublesome 
seasonal  dependence.  We  have,  therefore,  adopted  Pope's  basic  revision 
of  replacing  the  sunspot-number  dependence  in  our  old  high-latitude  term 
with  a dependence.  We  have  not  used  his  exact  formulation,  however, 
which  deleted  our  earlier  description  of  the  variation  in  statistical 
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width  of  the  scintillation  boundary  with  ita  location.  We  find  that  there 
is  such  a variable-width  behavior  in  the  in-situ  data  of  Sagalyn  et  al. 
(1974) ; it  is  described  by  the  denominator  of  the  error- function  argument 
in  Eq.  (44) . 

Having  incorporated  a K dependence  in  Eq.  (44),  it  was  a simple 

P 

extension  to  include  the  K^-de pendent  migration  of  the  auroral  oval  in 

Eq.  (45).  A more  accurate  index  to  use  might  be  the  auroral -electrojet 

index,  AE,  but  we  have  elected  to  use  only  for  simplicity.  The  two 

magnetic  Indices,  K and  AE,  are  rather  highly  correlated,  as  shown  in 

P 

Figure  4.  To  avoid  the  necessity  for  a user  to  specify  two  magnetic 
indices  rather  than  just  one,  we  have  elected  to  exclude  AE  from  the 
model  unless  later  analysis  of  scintillation  data  shows  that  including 
it  would  improve  the  models  accuracy  significantly.  In  a similar  vein, 
we  have  avoided  including  a geomagnetic  dependence  in  the  equatorial 
tens,  retaining  instead  the  original  sunspot-number  dependence.  Further 
experimental  work  is  necessary  before  a reliable  geomagnetic  dependence 
can  be  incorporated  in  the  equatorial  term;  a likely  candidate  for  a 
magnetic  index  is  D^,  the  seasonal  and  solar-cycle  behaviors  of  which 
are  similar  to  those  of  scintillation. 

In  addition  to  the  above  modifications  to  the  1973  model,  some  of  the 

constants  employed  have  been  changed.  These  changes  stem  partly  from 

improved  calibration  of  some  data  sets  used  in  development  of  the  original 

model.  The  data  sets  in  question  (Aarons  et  al.,  1963;  Koster,  1968; 

Joint  Satellite  Studies  Group,  1968)  were  presented  in  terms  of  a hand- 

scaled  scintillation  index,  SI  (Whitney  and  Malik,  1968) , which  only 

recently  (Whitney,  1974)  has  been  calibrated  reasonably  accurately  to  S^. 

It  is  still  to  be  noted  that  the  relation  between  SI  and  S,  must  be 

4 

scintillation-rate  dependent,  and  data  characterized  directly  by 
are  strongly  preferred  for  modeling. 
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AE,  AND  PLANETARY  MAGNETIC  INDEX, 


Figures  5 through  8 show  the  degree  of  fit  between  the  model 
now  in  RFMOD  and  several  data  sets  that  were  presented  by  the  observers 
either  in  terms  of  or  in  terms  of  an  empirical  index  such  as  SI, 
which  we  then  converted  to  S^.  Figures  9 and  10  compare  the  amplitude 
pdf  and  cdf  calculated  by  means  of  DIST  with  corresponding  histograms 
developed  from  ATS-5  data  provided  by  Aarons  (private  communication, 
1974).  The  corresponding  pdf  for  phase  is  shown  in  Figure  11;  no  phase 
data  exist  at  this  time  for  comparison. 

IV  THE  CHANNEL -MODEL  CODES 

Of  the  two  computer  programs,  RFMOD  and  DIST,  the  former  is  by 
far  the  more  complicated.  It  is  a modification  of  the  program  BPMOD, 
described  by  de  la  Beaujardiere  and  McNiel  (1971)  , and  consists  of  the 
routines  listed  in  Table  1.  Its  very  simplified  flow  chart  is  shown  in 
Figure  12.  Using  RFMOD  consists  essentially  of  preparing  input  cards, 
of  which  there  are  three  types:  (1)  a title  card,  (2)  namelist  cards, 

and  (3)  a changing-parameter(s)  card.  Three  sample  sets  of  input  cards 
are  shown  in  Figure  13,  and  each  type  of  card  is  described  below. 

The  title  card  is  used  to  assign  any  desired  title  of  up  to  80 
characters  to  a calculation  case  or  series  of  cases.  It  is  read  by 
the  driver  program,  RFMOD,  into  the  array  LBL  with  a 20A4  format,  and 
printed  in  the  output  by  the  subroutine  IMPRIME;  it  has  no  effect  what- 
soever on  the  rest  of  the  program.  Any  number  of  cases  can  be  run 

under  a common  title,  which  will  be  printed  in  the  output  for  each  case. 

The  namelist  cards  are  used  to  input  the  fixed  parameters  for  a 

calculation  case.  They  are  read  by  subroutine  READIN,  by  means  of  the 

namelist  INPAR.  The  user  gives  the  name  and  desired  value  for  each 
quantity  listed  above  the  dashed  line  in  Table  2 that  is  to  be  fixed. 

The  number  under  the  heading  K in  Table  2 is  an  index  used  in  the 
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INTENSITY-SCINTI LLATION  INDEX  OBSERVED  *00-MHz  INTENSITY-SCINTILI  ATION  INDEX, 

(Xs)  BY  PREDDEY,  MAWDSLEY,  AND  IRELAND  FOR  THREE  STATES  OF  GEOMAGNETIC  DIS- 

• 1869)  IN  THE  MID-LATITUDE  REGION  TURBANCE,  OBSERVED  (symbol*)  BY  EVANS 

(1973)  NEAR  THE  SUBAURORAL  SCINTILLATION 
BOUNDARY 


FIGURE  9 COMPARISON  OF  CALCULATED  (smooth  curve) 
PROBABILITY  DENSITY  FUNCTION  FOR 
AMPLITUDE  WITH  OBSERVATIONS  (histogram) 
OF  ATS-6  TRANSMISSION  AT  137  MHz.  MADE 
BY  AARONS  (private  communication,  1974) 
THROUGH  THE  MID-LATITUDE  IONOSPHERE 
FROM  NEAR  BOSTON 


FIGURE  10  COMPARISON  OF  CALCULATED  (smooth  curve) 
CUMULATIVE  PROBABILITY  FUNCTION  (cpf) 
FOR  INTENSITY  (or  for  amplitude)  WITH 
OBSERVATIONS  (histogram)  FOR  CASE  IN 
FIGURE  9.  (Calculations  of  cpf  art  intended 
for  use  in  assessing  fade-margin  requirements.) 
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FIGURE  11  CALCULATED  PROBABILITY  DENS'TY  FUNCTION  FOR  PHASE  FOR  THE 

CASE  IN  FIGURES  9 AND  10.  (No  observed  phase  data  exist  for  comparison.) 


Table  l 


ROUTINES  IN  PROGRAM  RFMOD 


Name 

Function 

BLOCK  DATA 

Initializes  common  parameters. 

SUBROUTINE  READ IN, 
with  ENTRY  INC 

Performs  program  inputting,  and  increments  the 
changing  parameter(s) . 

SUBROUTINE  SCINT 

Controls  the  basic  scintillation  calculation  and 
provides  quantities  needed  for  SUBROUTINE  SIGS. 

SUBROUTINE  SIGS 

Performs  the  basic  scintillation  calculations. 

FUNCTION  RMSDN 

Contains  the  worldwide  model  for  strength  of 
F-layer  irregularities  used  in  the  scintillation 
calculation. 

FUNCTION  RFBPT 

Evaluations  the  diffraction  integral  by  an 
algebraic  approximation. 

SUBROUTINE  FINDZ 

Computes  the  reduced  distance  from  scattering 
layer  to  receiver  for  use  in  evaluating  Fresnel- 
ao ne  radius. 

SUBROUTINE  MGFLD 

Contains  an  earth-centered , axially  tipped, 
dipole  model  of  the  geomagnetic  field  and  cal- 
culates the  geomagnetic  coordinates  of,  and  the 
field  components  at,  a specified  point. 

subroutine  aznside 

Given  the  latitude  and  longitude  of  two  points, 
A and  ft,  computes  the  azimuth  of  B as  seen  from 
A and  i.he  great-circle  arc  AB. 

SUBROUTINE  COORD 

Calculates  coordinates  of  a desired  point  on  a 
specified  great-circle  arc. 

FUNCTION  ERF 

Computes  the  Error  Function. 

SUBROUTINE  CFRMS 

Controls  calculation  of  rms  phase. 

FUNCTION  PHAS, 
with  ENTRY  PRT 

Evaluates  the  Hatfield  phase  distribution  of 
a scintillating  signal. 

I UNCTION  FI,  with 
ENTRY  F2 

Provides  integrand  functions  for  evaluating 
moments  of  phase  distribution. 

SUBROUTINE  SIMPA 

Adaptively  performs  numerical  Integration 
using  Simpson’s  rule. 

SUBROUTINE  IMPRIME 

Performs  program  outputting. 

LA-1078-28H 


FIGURE  12  SIMPLIFIED  FLOW  CHART  OF  PRO-  HAM  RFMOD 


34 


V 


I 

r: 


/ 

$ 

4i 

*• 

I 

> 

i 

! 

• 

i 

I 


o 

o l 

is 

Ui  UJ 

1 1 

H 2 


w 

E 

£ 

1 

E 

2 

w 

z 

5 _ 

H 

55 


« JO 


• I *- 

n|  n 

-N  ■ kn 

•»  3 1- 

«•  t w*8 

• ic  — 

NR  n 

0 fCirt 

•0  P 04 

n e «>c 

oie- 

N R n 

♦ R m 

r* 

■ Km* 

oR  - 

4»n  R n 

♦ R M* 

<0  R o* 

••  P •»  R 

IM*  « 

^ R kn 

«0  a o- 

mm  P 0i  k 

o * »- 

IN  IC  ♦*» 

N«tn 

40  P 04 

40  P 0*  p 

o * «— 

ON  * C.J 

•0-  7t  fcT» 

40  a 04 

0*  cnX 

4bi  a *— 

NC  »4 

Nr  R in 

40  e 04 

«•  n •*  p 

o C »- 

on  ti  o» 

ir  C in 

40  P 04 

•»  p to*  P! 

e C ^ 

M B « 

NT  K U- 

40  K 04 

« B 0C 

• IT  — 

IN  H 14 

n S 

«0  A r* 

0 P 0* 

ON  8 *0 

no  ft  if* 

40  * 04 

mm  8 408 

0 3- 

N J W 

NT  ft  U» 

40  a 04 

« 8 03 

« & — 

IN  S «*» 

n ; o 

40  S 04 

«0  & 40  & 

o » — 

N 3 14 

no  » vr* 

0 «|4 

Ml  0| 

m»  8 *-* 

fs  S 44 

N B in 

wfirs 

•0  a •*« 

«=>  X — 

Mj  14 

o 3 m 

40  a 04 

m3  03 

«»  8 — 

NOn 

•0  Q *« 

40  Q 04 

mm  S 0*0 

w 3 — 

cm  S n 

o ttw 

« S'** 

mS0O 

N 5 14 

o Km 

<0  6 14. 

m 5 MS 

Pl- 

IN.  S 14 

o 3 m 

40  |N 

••  a m s 

o3- 

4N  ft  44 

oSui 

« Sn 

m a m s 

«»  » — 

IN  X 44 

o S m 

40  a o- 

m a 0i 

o K •— 

IN  R 44 

N Sin 

40  U O* 

M R 0R 

• * — 

IN  31  »4 

o R in 

40  a o— 

Ml  0| 

IN  SI  14 

0*  8 if* 

o Bn 

mm  TX  40« 

• X — * 

IN  * *4 

oSo 

mji- 

M 3 0]t 

«»  3 ~ 

IN  R *4 

4 a in 

*•  a *4 

••  a 0*  a 

«•  3 *• 

IN»I4 

♦ »« 

40  a o— 

«•  a or 

IN  R 14 

o Km 

40  5*  o- 

0»  S M>  K 

•*  5 — 

IN  31  14 

N Jm 

0*0 

M*m*o 

IN  C 14 

4T  3 *n 

0*0 

-*«s5 

O « - 

IN  « 14 

o 9 m 

40  3 o- 

m ” 4 *u 

» V — 

N 5 44 

*0  C »* 

I#  5 r- 

— 5 —^o 

«D|- 

in  f n 

40  m 

40  3 04 

- * - "'•! 

«•  c — 

IN  # n 

■0“  « 0» 

40  V*  O- 

••  c 0*x?gj 

mm  * — ’ 

IN  » 14 

^ 3 m 

0 JO 

•0  3 •*  360 

• 3 — 

IN  44 

0*  3 a* 

09f 

•0  « 0*U* 

ON  S O'* 

0-  » *» 

080 

• 9»!5I 

«»;*•— 

IN  5 »4 

0-  « if* 

10  4 04 

«0  4 0»  « 

«»  3 •— 

IN  ? 14 

o3m 

0*0- 

— 3 0*rf 

«=>  R — 

IN  S *4 

<0  ft  •/* 

40  T,  o— 

•Rf*5, 

IN  ^ 14 

1 *m 

0*0- 

— a 0*  *5 

40  0 

IN  ^ 14 

o R m 

40  R 0>4 

•0  R 40  KU 

•»  g 

IN  «I4 

4T  ft  m 

40  * f 

••  5 •**4 

o R *“ 

ON  P 0* 

o P.  m 

0*1^ 

«0  « 0*  «g} 

o X ** 

IN  * **» 

o S m 

0*0- 

— X 0*Xp 

n — — 

N B 14 

o B m 

40  a 04 

M fl  0*  « J 

•»  « — 

on  » «•* 

o Bm 

40  a 04 

0B«fl2 

N R 14 

40  C IT* 

40  n 04 

0R0  Fi 

o x — - 

* ** 

Vlw 

0*o- 

•0  a 0*a 

«— » •— 

IN  R *4 

PR» 

0 R ^ 

M R 0R 

C*  R — » 

IN  H 14 

0 * U* 

40  a 04 

M H MR 

» fi 

on  r o» 

o am 

40  R 04. 

MM  R 0*  R 

IN  « 14 

40  a if* 

0*0* 

0*  a •»  a 

a UJ  «- - 

CN  R 14 

40  a m 

40  « 04 

« 0»  a 

«f  X — 

w R *4 

40  X If* 

40  a r4 

«•»•*« 

IN  B 44 

n R m 

40  n 04 

«0  & Mt  8 

o n — 

IN  R 44 

o a m 

0RN 

•0  n •*  n 

<0  R — 

444  R O* 

0 R 0* 

0RO. 

MM  R 40  R 

«»  R — 

IN  ft  *44 

40  * If* 

0*o- 

•m  a o*  a 

a 3 

IN  » 14 

■O  * in 

40  * 04 

MM  » 0*  * 

40  0 •“ 

IN  P 14 

■4-  * U4 

40  1*  04 

«f  * M* 

o P •— 

IN  B 14 

0 R if* 

40  R 04 

M B OtC 

O * — * 

IN  * 14 

40  0 If* 

40  * 04 

M R Ol  P 

o r - 

IN  « 14 

40  V If* 

40  IP  04 

M»  » 4 0 * 

«s  I 

IN  J 14 

O * m 

40  X 04 

M X 0* 

40  O *» 

IN  C O 

•O  Q if* 

40  e i4 

«M  a «n  o 

0 2 — 

r?  <•» 

o c m 

40  1— 

MO  t-  0» 

o r;  - 

ON  - «** 

N c m 

<0  SB  1— * 

•M  E 0>  1=  O 

C3  « r— 

!N  2 »4 

o 2 m 

<0  S r- 

05W5* 

ca  « *— 

IN  • 14 

40  0 If* 

U>  •»  04 

MM  0 40  0 ^ 

CJ  « — 

IN  “ 14 

o • in 

40  0 f- 

M • >n  «n 

ON  - 04 

40  - u* 

40—04 

mm  — o>  — 

O W — 

04  • 04 

«0  <•  IT 

40  «*  04 

0 N . M 

o « — 

ON  w.  14 

0 **  * 

40  tf*  |4 

M HI  0« 

ON  - 

o * in 

m n r- 

0 * W* 

_ #f,  n 

O n in 

40  r» 

M n 0 - 

a r,  - 

IN  4 14 

o 4 m 

— (4 

0 «■  0 N 

CT>  - — 

— C4 

40  ■—  If* 

U»  — 04 

«o  — «n  — 

I 


i 


i 

t 


i 


FIGURE  13  SAMPLE  INPUT  CARDS  FOR  PROGRAM  RFMOD 


Table  2 


DATA -CARD  VARIABLES 


K 

Variable  Name 

Unit 

Definition 

1 

FREQ 

Megahertz 

Frequency 

2 

SSN 

None  (floating  point) 

Sunspot  number 

3 

DAY 

None  (floating  point) 

Day  of  year 

4 

TIME 

Decimal  hours 

Meridian  time  at  iono- 
spheric penetration  point 

5 

RLAT 

Decimal  degrees 

Receiver  latitude 

6 

RLON 

Decimal  degrees 

Receiver  longitude 

7 

HR 

Meters 

Receiver  altitude 

8 

TLAT 

Decimal  degrees 

Transmitter  latitude 

9 

TLON 

Decimal  degrees 

Transmitter  longitude 

10 

HT 

Meters 

Transmitter  altitude 

11 

FKP 

Decimal  value  0.0 
through  9.0 

Planetary  magnetic 
index,  Kp 

-- 

11 

None  (value  l or  0) 

Data  flag  (see  text) 

12 

RCRD 

Decimal  degrees 

Receiver  coordinates 

13 

TCRD 

Decimal  degrees 

Transmitter  coordinates 

common  VARPAR  to  identify  the  parameters.  The  user  may  list  the  fixed 
parameters  in  any  order  on  the  namelist  cards. 

A parameter  input  on  a namelist  cart*  will  keep  the  specified  value 
until  changed  by  means  of  a new  namelist  card  or  until  input  on  a changing- 
parameter(s)  card.  The  first  column  of  each  namelist  card  must  be  blank, 
and  the  first  characters  on  the  first  namelist  card  must  be  &INPAR.  The 
last  namelist  card  in  a sequence  must  terminate  with  the  characters  &END. 
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The  changing-parameter(s)  card  is  employed  by  the  user  to  select  the 
independent  variable(s)  for  a given  calculation  case.  He  may  choose 
either  of  two  program  modes.  In  Mode  1,  he  chooses  any  one  of  the  param- 
eters listed  in  Table  2 for  which  K is  between  1 and  11,  inclusive,  and 
inputs  its  name  on  the  changing-parameter  card.  In  this  mode,  the  first- 

order  signal-statistical  parameters,  including  S.  and  $ , are  calculated 

*+  rms 

as  functions  of  the  single  changing  parameter.  In  Figure  13,  the  cards 
identified  as  Set  1 are  for  two  sequential  cases  run  in  Mode  1,  under  a 
common  title. 

Mode  2 permits  the  user  to  increment  the  position  of  either  the  re- 
ceiver or  the  transmitter  along  an  arc  of  a great  circle.  In  this  mode, 
both  the  latitude  and  longitude  of  the  receiver  or  of  the  transmitter  be- 
come changing  parameters,  and  the  signal-statistical  parameters  are  cal- 
culated for  equidistant  points  along  the  arc.  Mode  2 is  selected  by  naming 
RCRD  or  TCRD  on  the  changing-parameters  card  (see  Table  2) . In  Figure  13, 
the  cards  identified  as  Set  2 are  for  a run  in  this  mode. 

The  required  contents  of  the  changing-parameter(s)  card  are  listed 
in  Table  3.  Note  that  in  Mode  1,  if  the  final  value  of  a changing  param- 
eter is  smaller  than  the  initial  value,  the  increment  value  must  be  speci- 
fied as  negative.  Absence  of  the  required  negative  sign  in  this  situation 
will  cause  the  run  to  terminate,  and  an  error  message  will  be  printed. 

After  reading  a set  of  cards,  the  program  will  keep  the  same  title 
and  read  only  namelist  and  changing-parameter(s)  cards,  so  long  as  II  * 1. 
After  encountering  II  * 0,  it  will  read  a new  title  card.  The  program 
will  terminate  if  the  name  on  the  changing-parameter(s)  card  does  not 
correspond  to  one  of  those  numbered  1 through  13  in  Table  2.  Thus,  to 
terminate  a run,  the  user  may  prepare  a changing-parameter  card  contain- 
ing an  unfamiliar  word,  such  as  the  FINI  appearing  in  Figure  13. 


* 


37 


Table  3 


CONTENTS  OF  CHANGING -PARAMETER  S)  CARD 


Columns 

Mode-1  Content 

Mod-^-2  Content 

Symbol 
in  Program 

1 through  4 

Name 

Name 

Name 

9 through  16 

Initial  value 

Initial  latitude 

Al 

17  through  24 

Final  value 

Initial  longitude 

A2 

25  through  32 

Increment  value 

Final  latitude 

A3 

33  through  40 

-- 

Final  longitude 

A4 

41  through  42 

-- 

Number  of  increments 

NA5 

For  model  development  and  testing,  we  employed  a plotting  routine 
in  conjunction  with  RFMOD  to  produce  graphs  such  as  those  in  Fibres  5 
through  8.  The  available  version  of  RFMOD  provides  an  array  that  the  user 
may  employ  with  his  own  computer  peripherals . In  the  absence  of  user 
modification,  the  array  will  be  output  in  tabular  form  along  with  identify- 
ing information,  as  shown  for  a test  case  in  Figure  14,  In  the  case  sh<^wn, 
the  changing  parameter  was  the  latitude  of  the  transmitter,  which  was 
taken  to  be  in  a polar-orbiting  satellite  passing  over  a receiving  station 
near  Boston  (namely,  the  Sagamore  Hill  site  of  the  AFCRL  group  headed  by 
Aarons) . The  changing  parameter  is  followed  in  each  row  by  the  azimuth 
and  elevation  of  the  satellite,  the  three  probability-ellipse  parameters, 
and  then  by  the  calculated  scintillation  indices.  The  input  cards  for 
this  run  are  identified  as  Set  3 in  Figure  13. 

It  is  envisioned  that  a user  will  inspect  the  values  of  S4  and/or 
PRMS  predicted  for  his  communication  link  by  RFMOD  and  then  decide  if  he 
wants  more  detailed  signal-statistical  information  for  a particular 
situation.  In  the  example  of  Figure  14,  for  instance,  he  might  decide 
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that  he  mats  the  full  first-order  distributions  of  intensity  and  phase 
for  the  case  in  which  S4  on  his  channel  is  .538  and  PBMS  is  32.835  degrees. 
To  obtain  the  distributions  he  would  employ  DIST,  inputting  RZRO  “ .371, 
BMAGR  - .924,  and  DELTA  - 3.419  degrees,  on  a single  card  in  F10.3  format. 
The  DIST  code  consists  of  the  routines  listed  in  Table  4. 


Table  4 

ROUTINES  IN  PROGRAM  DIST 


Name 

Function 

SUBROUTINE  PDFI 

Controls  calculation  and  output  of  probability 
density  function  for  amplitude,  by  means  of  one- 
dimensional integration  over  two-dimensional 
Gaussian  probability  density  function. 

SUBROUTINE  PARRAY 

Initializes  amplitude  array  for  SUBROUTINE  PDFI. 

SUBROUTINE  GAUP 

Controls  calculation  of  Gaussian  probability 
density  function. 

FUNCTION  PR 

Computes  two-dimensional  Gaussian  probability 
density  function. 

FUNCTION  XINTGR 

Performs  integration  by  Simpson's  rule. 

SUBROUTINE  PDFP 

Controls  computation  and  output  of  probability 
density  function  for  phase,  calling  SUBROUTINE 
PHAS. 

SUBROUTINE  PHAS 

Evaluates  Hatfield  phase  distribution. 

FUNCTION  ERF 

Computes  the  Error  Function. 

SUBROUTINE  CPFI 

Calculates  and  outputs  the  cumulative  probability 
density  function  for  amplitude  or  intensity 
(expressed  in  dB)  by  performing  Simpson's  rule 
integration  over  probability  density  function 
for  amplitude. 

Again,  for  modeling  we  employ  a plotting  routine  in  conjunction  with 
D1ST  to  produce  graphs  such  as  those  in  9,  10,  and  11.  The  available 
version  of  DIST  simply  outputs  the  array  from  which  such  graphs  are  made, 
in  the  line-printer  formats  shown  in  Figures  IS,  16,  and  17.  As  with 
RFMOD,  a user  may  wish  to  employ  DIST  in  conjunction  with  his  awn  plotting 
routine  or  in  some  other  output  form.  Alternatively,  he  may  wish  to 
reformat  the  array  listing,  perhaps  outputting  fewer  values  and  possibly 
presenting  CPF  in  terms  of  percent  or  decimal  link  reliability  (the  one's 
complement  of  CPF)  rather  than  directly  as  in  Figure  16. 

V CONCLUSION 

It  has  been  our  attempt  to  develop  a channel  model  that  may  be  used 
conveniently  for  calculating  first-order  signal  statistics,  expressed  in 
a format  relevant  to  systems  application.  The  resulting  codes,  RFMOD 
and  DIST,  represent  the  first-year  output  of  a broader  effort  to 
definitively  model  the  transionospherlc  communication  channel.  These 
codes  are  to  be  augmented  by  routines  to  permit  calculation  of  second- 
order  statistics,  such  as  fluctuation  spectra  of  amplitude  and  phase,  as 
well  as  signal  correlation  at  spaced  receivers  and  across  bands  of 
frequency.  Moreover,  the  existing  version  of  RFMOD  is  likely  to  evolve 
over  the  next  year  in  two  important  ways. 

First,  it  is  intended  to  modify  SUBROUTINE  SIGS  in  RFMOD  to  fully 
account  for  multiple  scatter.  Theoretical  work  now  underway  is  hinting 
at  a possible  form  for  the  modification.  The  single-scatter  theory 
described  in  Section  II  of  this  report  describes  the  simultaneous  rota- 
tion (increasing  5)  and  circularization  (decreasing  IBqI/Rq)  of  the 
first-order  probability  ellipse  for  the  received  signal's  complex  ampli- 
tude (illustrated  in  Figure  1),  as  the  wave  propagates  away  from  the 
scattered  medium.  In  the  single-scatter  limit,  the  ellipse  emerges  from 
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FIGURE  15  EXAMPLE  OF  AMPLITUDE  pdf  OUTPUT  FROM  AVAILABLE  VERSION  OF  PROGRAM  DIST 
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the  scattering  layer  highly  elongated  ( |bq|/R0  m 1)  and  oriented  nearly 
vertically  (6  » 0)  on  the  complex  plane. 

Multiple-scatter  calculations  for  the  special  case  of  normal  inci- 
dence on  the  scattering  layer  show  that  some  circularization  (|B0|/R0 
* « “*<>),  but  no  rotation  (&  « 0),  is  produced  by  the  rescatterings.  This 
is  an  intuitively  understandable  result  and,  if  corroborated  by  similar 
calculations  for  off-normal  incidence,  is  expected  to  lead  to  a multiple- 
scatter  correction  for  |B0|/Ro  analogous  to  Eq.  (28)  for  R^,.  It  will 
presumably  also  lead  to  a modification  of  Eq.  (24)  for  relating  the 
ellipse  orientation  to  its  elongation  at  the  receiver.  The  result  is 
certain  to  be  a reduction  in  the  third  term  of  Eq.  (5)  for  intensity- 

scintillation  index  as  the  fraction  of  flux  scattered  (R  ) increases. 

o 


It  will  presumably  limit  values  of  to  something  near  unity  and  improve 
agreement  between  calculated  and  observed  values  for  moderate  to  strong 
scintillation. 

The  other  Important  area  in  which  RFMOD  is  expected  to  evolve  is 
that  of  scintillation  morphology,  which  will  entail  updating  of  FUNCTION 
RMSDN.  As  more  definitive  scintillation  data  become  available,  the  model 
of  electron-density  irregularity  strength  contained  in  RMSDN  presumably 
will  be  improved  by  a process  that  may  be  viewed  as  calibration  [i.e., 
the  main  changes  should  be  adjustment  of  the  coefficients  in  Eqs.  (42) 
through  (45)].  This  will  be  particularly  so  in  regions  where  the  model 
is  currently  calibrated  in  terms  of  ad  hoc  indices  such  as  SI  instead  of 
in  terms  of  S^.  One  such  region  is  the  important  one  centered  on  the 
geomagnetic  equator.  Useful  equatorial  data  are  expected  from  observa- 
tions of  the  NQAA  ionospheric  beacon  on  board  NASA's  ATS-6,  particularly 
measurements  of  S^  at  140  and  360  MHz.  During  periods  of  weak  to 

moderate  scintillation,  estimates  should  also  be  obtained  of  $ 

7 rms 

at  140  MHz.  The  ATS-6  observations  also  will  provide  histograms 
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1 

of  amplitude  (under  all  conditions))  and  of  phase  (under  conditions  of 
weak  to  moderate  scintillation)  for  testing  the  predictions  made  by  the 
D1ST  code. 

Pending  the  improvements  described  above,  we  offer  RFMOD  and  DIST 
as  the  best  way  to  estimate  the  first-order  signal  statistics  to  be 
expected  at  the  output  of  a user's  transionospheric  communication  link, 
short  of  performing  a series  of  observations  over  that  particular  link. 

It  must  be  remembered  that  the  estimated  values  are  those  expected  for 
average  ionospheric  conditions  in  the  observing  situation  (frequency, 
transmitter  and  receiver  location,  time  of  day,  K , etc.)  specified  by 
the  user.  There  is  presumably  some  statistical  distribution  of  each 
parameter  estimated  about  its  calculated  value,  and  the  present  model 
does  not  account  for  such  statistical  variation.  Possibly  of  equal 
importance,  the  model  has  yet  to  be  subjected  to  testing  against  observa- 
tions above  400  MHz  (e.g.,  in  the  SHF  band). 
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